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We consider a three-dimensional weakly interacting Bose gas in the fluctuation region (and its 
vicinity) of the normal-superfluid phase transition point. We establish relations between basic 
. thermodynamic functions: density, n(T,/i), superfluid density n a {T,n), and condensate density, 

n cn d(r, fi)- Being universal for all weakly interacting |i/>| 4 systems, these relations are obtained from 
Monte Carlo simulations of the classical \i/j\ 4 model on a lattice. Comparing with the mean-field 
results yields a quantitative estimate of the fluctuation region size. Away from the fluctuation region, 
on the superfluid side, all the data perfectly agree with the predictions of the quasicondensate mean 
' field theory. — This demonstrates that the only effect of the leading above-the-mean-field corrections 

in the condensate based treatments is to replace the condensate density with the quasicondensate 
one in all local thermodynamic relations. Surprisingly, we find that a significant fraction of the 
density profile of a loosely trapped atomic gas might correspond to the fluctuation region. 



o 
o 

or 

CZ2 



o 

1 



- I - < 

>< 
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I. INTRODUCTION 



Typically, only special properties of critical systems (such as critical exponents and amplitude ratios) are universal 
(see, e.g. Ref. OJ. Weakly interacting U(l) models, including dilute Bose gases (BG), are an exception from this rule 
in the sense that all their thermodynamic properties are universal close to the critical point. This unique feature 
allows complete characterization of the fluctuation region and its vicinity for all such models without adjustable 
parameters. In the recent paper— devoted to two dimensional systems we obtained the universal thermodynamic 
' relations for the density, n, superfluid density, n s , and quasicondensate density, n qc , by performing high- precision 
l— ', Monte Carlo simulations of the classical |?/>| 4 model on a lattice. Here we present the results of a similar study of the 
| ■ three-dimensional (3D) case. 

Away from the critical point, the mean-field (MF) treatment of the weakly interacting BG^ is a reliable theoretical 
■ scheme controlled by the small dimensionless parameter 

n 1/3 a<l, (1.1) 



where a is the scattering length. In the pseudo-potential approach, one models the same scattering length a by 
introducing a weak short-range potential, U ps d, of radius r -C n -1 / 3 and amplitude U = J ^7 ps d(r) dr. The relation 
between U and a is provided by the perturbation theory (we set % = 1) 



U-^(l + tm f * i + (1.2) 



m \ J (27r) 3 p 



' 

The divergence of the second term in the ultra-violet limit is cut at p ~ 1/r; it cancels out when final answers are 
expressed in terms of the scattering length 3 -' 1 ' 5 . Close to the Bose condensation point, when n 2 / 3 ~ mT, one can 
write Eq. as 



In the fluctuation region, 
03 



mU^l/VmT. (1.3) 



|f| = |(T C - T)/T c \ ~ m 3/2 T 1/2 U , (1.4) 



(T c is the critical temperature), perturbative in U approaches do not WO rk5i£i£iSi2iifl. The physics now is determined 
by non-perturbative coupling between the long-wavelength components of the order parameter field, resulting in a 
crossover from MF to the generic U(l)-universality class critical behavior on approaching the critical point. 

An important observation about the fluctuation region of a weakly interacting BG is that in the small interaction 
limit all l^j 4 models — quantum or classical, continuous or discrete — allow a universal descriptio n 2 ^ 6 . The crucial 
circumstance is that long- wavelength components of the order parameter field, ij)(r), with momenta much smaller 
then the thermal momentum kx = \JmT are classical in nature (large occupation numbers), and only components 
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with k < fc c = m 2 TU <C &t are strongly coupled. In this limit the effective Hamiltonian is given by the classical \i/j\ 4 
model 

m = J ^m 2 + jM 4 -M\ 2 }dr, (1.5) 

where /x is the renormalized chemical potential for long-wave modes. The microscopic physics of different models is 
important only at momenta k 3> k c , where the system behavior is perturbative and thus may be easily accounted for 
analytically. 

Our approach to the problem starts from observation that n s , n cn( j, and n — n c , are universal functions of the 
shifted chemical potential, fi — fi c , where fi c (T), n c (T) are the (model specific) critical point parameters. To find these 
functions one has to solve accurately any of the weakly interacting \ip\ 4 models, and we resort to Monte Carlo (MC) 
simulations of the classical |^| 4 model on a lattice using very efficient Worm-algorithm not suffering from the critical 
slowing downii. In 2D, the same idea was successfully implemented in Ref. |2|. Before that, quantum to classical 
mapping was used in Refs. M llflrj to determine the critical point parameters (both in 2D and 3D). 

This approximation-free numerical solution yields an accurate description of the system thermodynamics. In the 
critical region, characterized by known exponents of the U(l) universality class, our goal is to find amplitudes of the 
power-law dependencies for all quantities in question. By comparing to the results of the MF theory, we establish 
quantitative limits on the size of the fluctuation region. Away from the fluctuation region our accurate data unam- 
biguously distinguish between MF theories based on the notion of condensate and quasicondensate — the amplitude 
of the order parameter field at intermediate distances — in favor of the latter. 

In the past, most non-pertur bative calcula tions using renormalization group (RG) and l/TV-expansions (for the most 
recent work see, e.g., Refs. I7ll8lll3ll4ll5ll6|) concentrated on the critical temperature shift, AT C , and the scattering 
in the results was quite substantial. In the absence of small parameters controlling the accuracy of the answer, the 
knowledge of final results is crucial in discriminating between competing approaches and in developing better schemes. 
However AT C , or An c , is only one of many universal properties of weakly interacting systems. The critical chemical 
potential shift as well as superfluid and condensate density behavior in the critical regime are also universal. A reliable 
theoretical approach should be able to reproduce all of them, and our results provide corresponding benchmarks. 

The paper is organized as follows. In Sec. [H] we use the analysis of dimensions to cast thermodynamic relations for 
the weakly interacting gas in the universal form. Special attention is paid to the ultra-violet and infra-red procedures 
of the chemical potential renormalization. In Sec. Illll we render MF theory based on the quasicondensate density. In 
Sec. lIVI we introduce the classical \ip\ 4 model on a lattice and the simulation algorithm. Our results are presented and 
compared to the critical and MF behavior in Secs.|3 In Sec. IVII we discuss our results in the context of experiments 
with ultracold gases and make comparisons with some analytical approaches to the fluctuation region. 



II. UNIVERSAL RELATIONS FOR WEAKLY INTERACTING \ip\ 4 MODELS 

In the U — > limit one can present simple arguments for the typical energy and density scales responsible for the 
non-perturbative behavior at the critical pointSt^^. To find the momentum separating weakly and strongly coupled 
modes, k c , one considers the three terms in the Hamiltonian (|1.5|) and determines when all of them are of the same 
order of magnitude for modes k < k c , assuming that modes with k > k c are already taken into account in renormalized 
values of the Hamiltonian parameters. This leads to the estimates 



where 



k 2 Jm ~ |0| ~ nU , (2.1) 

is the long wavelength contribution to the total density, and jl is the effective chemical potential for low-energy modes. 
In 3D one has n ~ fc^n^, and since k c is separating strongly coupled long- wave harmonics from slightly perturbed 
short-wave ones, the order-of-magnitude estimate for n^ c may be obtained from the ideal system formula: 

T T 
kc k 2 /2m-fl ' 

Substituting this back to Eqs. (|2. 111 - 12. 2[) yields 



k c = m 2 TU , n ~ m 3 T 2 U , |0| ~ m 3 T 2 U 2 . (2.4) 
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FIG. 1: The second order diagram for the self-energy E(w = 0,p = 0). 



The critical values of the density, or chemical potential, are, of course, model specific. [We find it convenient to 
work in the grand canonical ensemble and keep temperature fixed; in Sec. El we explain how all the results are easily 
converted into more familiar temperature plots.] In fact, in 3D the major contribution to n c , and fi c , is coming from 
high momenta k 3> k c . However, this model specific contribution to density does not depend on interactions in linear 
order in U, and thus can be easily calculated analytically. Thus if we count density and chemical potential from their 
critical values, we expect the dependence of n — n c on fi — fi c to be universal because of its long-wavelength nature. 
In view of Eq. 1)2. 4)1 we write the equation of state in the dimensionless form as 

n - n c = (m 3 T 2 U) X(X) , (2.5) 

where X is the universal control parameter 

x - ■ < 2 - 6 > 

with the typical variation across the fluctuation region of order unity. [For the gas in a slowly varying external 
potential, V(r), the function X(X) describes the spatial density profile, provided fj, in Eq. i)2.6)) is replaced with 
/i — V(r).] Similarly, 

n cnd = (m 3 T 2 U) f (X) , (2.7) 

n s = (m 3 T 2 U)f s (X) , (2.8) 

Establishing universal — for all weakly-interacting l^l 4 models — functions A, /o, and / s , and explaining how they work 
for the quantum Bose gas is the main goal of this paper. 

First, we note that n c and /i c for the Bose gas have been already determined in previous MC studies2*i2iii. The 
interaction induced shift of the critical density is universal 

n c = - Cm 3 T 2 U , C = 0.0142(4) , (2.9) 

where is the critical density of the non-interacting system; n' ' = (mT/3.313) 3 / 2 for the Bose gas. [We quote the 
mean of the results presented in Refs. 9 and 10 which overlap within the error bars.] 

Before quoting the result for /u c , we would like to discuss its subtle structurei2*iii£. Linear in U terms are accounted 
for by the MF expression 2nU. There are two distinctive contributions to [i c which are quadratic in U . One is universal 
and comes from strongly coupled critical fluctuations. The other one is perturbative and comes from the second-order 
diagram (see Fig.^l for the self-energy, £( 2 )(w = 0,p — 0): 

E P) = -2U 2 T 2 [ f dkldk2 Y 



( 27r ) 6 7^ 2 \ iu) ^ ~ e fcJN S2 - £fc 2 ][«(w Sl + u S2 ) - e kl+k2 ] 
2 ' 1 1 ' /k, ' /kj [N kl N k2 -N kl+k2 (N kl +N k2 +1)] (quantum), (2.10) 



(2tt) 6 J J e kl + e k2 - e kl+k2 

where e k is the dispersion law, uj s = 27rTs is the Matsubara frequency, and N k = \e tk l T — 1] _1 is the Bose distribution 
function. The corresponding classical expression is obtained by considering only the zero-frequency term [s\ = S2 = 0): 

l?] 2U 2 T 2 f f dkidk 2 . , . ,, , 
S (2) = - Nfi / / — - — classical , (2.11) 
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The self-energy integrals are divergent. In the infrared limit the divergence is logarithmic 

fe «fc T m 3 T 2 U 2 



dk/k . (2.12) 
For the Bose gas (e^ = k 2 /2m) it is also divergent in the ultra-violet limit as a power law 

- -7^ / / rfkldk27Vkl+k2 « / * = -2nU Una f * ^ = E« . (2.13) 

(2tt) 6 7 7 e fcl + e fc2 - e kl+k2 (2tt) 6 7 Vp 2 V J {^) 3 p 2 ) uv V 7 

We cast the last expression in the form which immediately relates it to the definition of the pseudo-potential, Eq. (|1.2fl . 
i.e. in the sum 2nU + T,^ the ultra-violet divergences cancel each other for terms cx U 2 . Identically, the same result 
is obtained if in all final answers we simply substitute U — > Ana/m and use 

s (2) _ E (2) _ E W ( quantum ) , (2.14) 

as the second-order self-energy for the quantum Bose gas. In what follows we employ this well established trick 3 4 5 . 
Obviously, in classical models quantum corrections to U are absent. Correspondingly, there are no divergencies in 
£( 2 ) apart from the logarithmic one, and the original expression for the self-energy should be used: = £< 2 ). 

The logarithmic divergence l|2.12|l can be simply truncated at some infra-red energy scale, or regularizediSiit 
Whatever the procedure, as long as it is the same for all weakly-interacting models, the difference between the true 
value of [i c and 2n c U + is coming from long-wave modes and thus is universal, while the 5> 2 ) term is model 
specific. In this paper we introduce an explicit infra-red cutoff by adding a constant k — k 2 /2m — m 3 T 2 U 2 /2 to the 
dispersion law, i.e. ^ + k everywhere in Eqs. (|2.1(J|) and (|2.11|) . Then 

m 3 T 2 U 2 ( r- Q s\ m 3 T 2 U 2 ( k T \ , , 

S (2) = j In ( 5.9030 Vm 3 T[/ 2 ) = - 2 In 5.9030 -f- ) , (2.15) 

The numerical value of [i c has been determined in Refs. fTojlTt For the quantum gas it reads 

fi c = 2n ( ^U+ Hfl^ll in ( .4213(6) = 2n c U + S< 2 > - (m 3 T 2 U 2 ) 9 , (2.16) 

7T \ k c / 

where 

6 Q = (1/tt 2 ) ln(5.9030/0.4213) - 2C = 0.239(1) . (2.17) 

In the MF theory, the definition of the effective chemical potential jl involves subtraction of the self-energy contri- 
butions 

p, = fj, ~ 2nU - E (2) . (2.18) 

This quantity is model-independent, and we find it convenient to introduce the corresponding universal function 9(X) 

2nU + ± {2) - n = m 3 T 2 U 2 e{X) . (2.19) 

It is, of course, directly related to the equation of state function \(X) 

9(X) = 2X(X) - X + 9 . (2.20) 

If we were to change the infra-red cutoff in the definition of self-energy, k — > k' we would have to modify the value 
of 9 accordingly 9 -> + (1/27T 2 ) In(/c//c'). 



III. MEAN-FIELD THEORY 



Away from the fluctuation region our simulation is supposed to reproduce the known perturbative results for a 
weakly interacting Bose gas. That is we need 9, A, /o, and f s as functions of X, the dimensionless variable 1/X 
playing the role of the small parameter that guarantees the applicability of perturbative treatment. The leading 
terms in the perturbative expansion away from the fluctuation region are ~ \X\, the next-order terms are ~ 
We confine ourselves to considering only these terms, ignoring the higher order corrections (in fact, we are not aware 
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of existing results for them). Hence, our numerics is expected to agree with the known analytic results only at large 
enough \X\, and only up to some constants. 

This degree of accuracy can be achieved within the MF treatment (plus an extra care of the higher-energy logarithmic 
correction to the chemical potential) both on the normal and superfluid sides. In the superfluid region, however, it 
is essential, that MF be based on the quasicondensate, rather than genuine condensate. In the condensate-based 
treatments the same accuracy is achieved only if beyond-thc-MF corrections are taken into account. The logarithmic 
correction to the chemical potential, & 2 \ has been already discussed in Sec.|H] For practical purposes, this correction 
is not large — unless the gas parameter n^^a is exponentially small, and taking it into account really makes sense only 
if all the ^-independent terms are accounted for as well. Theoretically, however, this correction involves the ultraviolet 
physics and thus is model-specific. Hence, if we ignore it, then we cannot render our answers model-independent. On 
the good side, we do not need to calculate S^ 2 ^ more accurately than it was done in Sec. [HJ since that would lead 
to the higher-order corrections which we ignore here. Thus, we simply add the term (|2.15|) to MF expression for the 
chemical potential of the Bose gas. 



A. Asymptotic behavior in the normal phase (X 



-oo) 



In the normal region X < we have /o and f s identically zero. Hence, the only quantity we have to look at is 
the equation of state X(X) or 9(X). Using MF equation for the effective chemical potential fi = (J, — 2nU ~ S'- 2 -' = 
—Om 3 T 2 U 2 , we calculate the difference between the density of the interacting gas n(T, jj) and the critical density of 
the ideal gas keeping only the leading linear in U terms: 



dk 



T 



9m 3 T 2 U 2 



Vdm 3 T 2 U 
V2tt 



We now notice that n 
derivation 



(0) 



[C - X(X)]m 3 T 2 U, see Eqs. £3 and O, and use Eq. 



V2 



— Ve = 2C + 6(0) -X at X 



(3.1) 

to complete the 
(3.2) 



B. Mean-field description of the superfluid region (X — > oo). Quasicondensate 

The standard MF approach to the weakly interacting Bose gas deals with the condensate density, rt cn d, defined 
through the one-particle density matrix, 

p(r) = (^(r)V(O)) (3.3) 

[ip is cither the field operator (in the quantum system) or a complex valued field (in the classical system) ; in the latter 
case tpi = tp*], as 

ricnd = lim p(r) . (3.4) 

r — >oo 

Well inside the region of applicability of the MF description, but not too far from the fluctuation region, the MF 
theory based on n cn a is less accurate than the theory dealing with the notion of gwasi-condensate. Of course one 
can go beyond the MF description in the condensate based theory and evaluate the corresponding corrections. It is 
important, however, that the quasicondensate MF description automatically captures the first order corrections to 
the condensate MF. [In 2D the quasicondensate MF has demonstrated perfect agreement with MC simulations away 
from the fluctuation region^.] We thus find it instructive to resort here to the quasicondensate MF description. 

The notion of quasicondensate was introduced by Popov^ ( "bare condensate" in the original version) . Basically, it is 
used to describe the order parameter with fluctuating phase (see alsoi^) and implies the possibility of parameterizing 
the field ^(r) for the weakly interacting system as 

V(r) = Vo(r)+^i(r) , (3.5) 
!fo(r) w ^V e ** (r) > ( 3 - 6 ) 
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where n qc is the quasicondensate density, and ipi is the Gaussian field independent of i^q. The Gaussian field ip\ is 
primarily responsible for the decay of p(r) at distances comparable to the thermal wavelength. The relation between 
n cnc i and n qc immediately follows from (|3.5(1 - H3.6|1 : 



^cnd 



n qc lim ( e »W-»(o)) 

r — *oo 



(3.7) 



When fluctuations of the phase field, ^(r), become noticeable — this is precisely the case in the vicinity of the critical 
point — the difference between n cn( j and n qc should be taken into consideration as well. 

There are strong arguments that it is n qci and not n cn( j, that is relevant to all the characteristics of the system, 
with n cn d being just one of them. Indeed, the physics at large distances, or small, but finite momenta, including 
the spectrum of elementary excitations, is determined by what is happening on smaller lengthscales, in other words, 
"high-energy" physics ultimately determines what happens at lower energies. The condensate is the macroscopic 
characteristic of the system; it should be derived from other finite-A: properties, and n qc governs them. 

This fact has been rigorously shown by one of us2S. Though the analysis of Ref. |2(] has been done for 2D systems 
(where the notion of quasicondensate is of crucial importance and cannot be avoided), it is applicable to the 3D case 
as well. The actual treatment closely follows Popov's theory^, with the only (but important) exception that n qc is 
understood and treated as a physical quantity. [Popov treated bare condensate as an auxiliary mathematical quantity 
explicitly dependent on the momentum k' separating 'slow' harmonics from the 'fast' ones. Correspondingly, he 
attempted to exclude this quantity from final answers to render them fc'-independent. — At low enough temperatures 
this can be easily done by replacing n qc — > p,/U . However, in the high temperature region the quantity n qc (T) is a 
meaningful physical parametefli^.] 

One may consider Eq. (|3.6|) as the definition n qc — it is the modulus of the order parameter field at large distances. 
This quantity appears in all MF equations just like n cn d does at low temperature. Then, n qc and T can be chosen as 
convenient independent thermodynamic parameters specifying the state of the system, the rest of the characteristics 
being expressed as functions of (n qc ,T). The basic results are as follow a 19 i 20 : 

n = n qc + n' , (3-8) 



with the non-quasicondensate part of the particle density, n', given by the integral 



n = 



d d k 
{2n) d 



where e(fc) = k 2 /2m is the free-particle dispersion 



e(fc) + n qc U - E(k) e(k)N E 
2E(k) E{k) 

i law, and 



E(k) = ^e{k)[e{k) + 2n qc U] 



(3.9) 



(3.10) 



is the Bogoliubov quasiparticle spectrum. 

For the one-particle density matrix one obtains 2 ^ (see also Ref. for a numeric check in 2D) 



p(r) 



p{r) 



(3.11) 



A(r) = 



d d k UN E 
— [l-cos(k.r)] 



E 



(3.12) 



p{r) 



^qc 



d d k 
(2-K) d 



cos(k-r) 



-n qc U- 



2E 



E e Ne 



E 



(3.13) 



C/ 3 / 2 after integration, and may be omitted 



[In Eq. H3.13fl , as well as in Eq. (|3.9|) , the first term in square brackets is 
in the region T 3> nU addressed in the present paper.] 

That n qc is the relevant quantity behind the long- wave physics is clear from the structure of Eqs. H3.11fl - H3.13fl . 
The second term in Eq. (|3.13l) decays first as a power law ~ (^r)" 1 , and then exponentially fast, so that on large 
length-scales the amplitude of the order parameter is given by n qc . Phase fluctuations are described by A(r). They are 
negligible at short distances of order l/fcy, but at distances ^ 1/ ^/mn qc U their contribution to the density matrix 
results in the difference between n qc and n cn( j, see Eq. (|3.7|l : 



^-cnd 



-A(oo) , 



(3.14) 
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The trick to solve MF equations in the vicinity of the critical point (1<I< fcr/^c) is to calculate differences 
between the corresponding densities (this makes all integrals to converge at energies <C T) keeping only the universal 
leading low-momenta terms. To get n' we add and subtract the ideal gas critical density 

rn*T 2 UJJZ 



n<°> - n' 



f d 3 k 


T eT~ 


1 (2^) 3 


7 ~ E 2 ~_ 



where 



> 1 (at X > 1) 



(3.15) 
(3.16) 



/qc m 3 T 2 U 

is the dimensionless parameter of asymptotic expansion for all thermodynamic quantities away from the fluctuation 
region. 

Combining l|3.15|l with (|3.8|) yields 



n = n^> + m T U [f qc ~ y/f qc /n] , (3.17) 
and with i|2.9[) we get 

A= / qc - y/U/n + C at X^oo. (3.18) 
The relation for the condensate density immediately follows from the asymptotic value of the phase correlator 

d 3 k UT 1 



A(oo) 



(2tt)3 E\k) 2n Jf q 



Hence 



fo=Ue~^VU) at X 



30 



To find f s , we consider the standard expression for the normal component density 4 



1 



d 3 fc 



3m J (2tt) 3 



dN E 



dE 



Integrating by parts, we rewrite Eq. (|3.21[1 as 

_ 2(2m) 3 / 2 
,ln ~ 3(2tt) 2 

We now subtract from n n the non-condensate density to evaluate the integral explicitly, and use identity n n — n' 
n-qc - n s to get 



u 



(3.19) 
(3.20) 
(3.21) 

(3.22) 



(3.23) 
(3.24) 



fa = fqc - \//qc/ 37r at X —> 00 • 

Finally, the quasicondensate density can be related to the chemical potential as 1 ^ 

H = {n qc + 2n')U + £ (2) = 2nU - n qc U + £ (2) , 

where we have added the term in accordance with the previous discussion. In contrast to the analogous relation 
in terms of the condensate and non-condensate densities^, this relation does not imply corrections of the order \f~X. 
Comparing Eq. i|3.24[l with the definition of 8- function, Eq. I|2.19[l . we see that 

9 = / qc at I^oo, (3.25) 



qc 



and with the 9-X relation H2.20JI and Eq. Ij3.18|l for A we get a self-consistent equation for /, 

oo . 



/ qc - 2^/qc/TT = X - 9(0) - 2C at X — ► oo . (3.26) 

Equation IpOIIjl along with Eqs. <|3~T%|) . ftTHfy . and (j3~2~3")> define 9, A, / , and f s as functions of X. 

We are in a position now to demonstrate that quasicondensate MF reproduces results of the condensate-based 
diagrammatic technique which accounts for leading, ~ y/X, corrections to the condensate MF answers^. Indeed, in 
the limit X — ► oo, Eqs. I|3.26|l and (|3.20(l define the effective chemical potential dependence on the condensate density 
as 



X w fo(l + Vfo/2w) - 2sJf /Tr = f - ixjhl2ix . 

Similarly, 

fs ~ /o(l + V^/2vr) - ^/37r = f + y/h/en , 
in complete agreement with Ref. 5. 



(3.27) 
(3.28) 
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FIG. 2: The #(X)-function and the isothermal density profile X(X) along with the corresponding MF expressions, Eqs. 
and 13. 251 . shown by dashed lines. 
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0.25 



A 



0.2 



X 




FIG. 3: The condensate and superfluid density dependencies on X and the corresponding MF predictions based on Eqs. I|3.2()|l 
and 13.23H shown by dashed lines. The initial parts of the plots are shown in the insets. 



IV. LATTICE ji^l 4 MODEL AND SIMULATION ALGORITHM 

As explained in Sees. Ill and ITT1 one can introduce universal functions for any model described in the long-wave limit 
by the Hamiltonian 1|1.5|) with small U. Classical lattice models are easier to deal with numerically, and there are very 
efficient classical algorithms which allow simulations of very large system sizes. In the present study, we performed 
simulations of systems with up to 128 3 lattice points. In addition, results obtained for a classical model directly test 
the idea of universality, since they have to agree with all MF predictions formally derived for the quantum Bose gas. 

Our simulations were done for the simple cubic lattice Hamiltonian 



h= Yl [£(k)-M]h/> fc | 2 + ^£h/#' 



(4.1) 
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where ipk is the Fourier transform of the complex lattice field ipi, and 

£k = —z V [2 - cos(fc Q a)] , (4.2) 
ma z z — ' 

a—x^y^z 

is the tight-binding dispersion law with momentum k confined within the first Brillouin zone (BZ). The self-energy 
expression H2.11JI was evaluated with this dispersion law and k = k c /2m numerically for each system size, and used 
then in the definition of the function 6(X), Eq. (|2.19|) . 

We employed the Worm algorithm for classical statistical systems^ which has Monte Carlo estimators for all 
quantities of interest here and does not suffer from critical slowing down. In cluster methods^ (which have comparable 
efficiency) the calculation of the superfluid density is more elaborate. 

From simulations of a series of system sizes and coupling constants, U = 2,1, 0.5, 0.25, we eliminated finite-size and 
finite-?/ corrections to the final results for 9(X), X(X), fo(X), and f s (X). 



V. SIMULATION RESULTS 



In Figs. [3 and |21 we present the final outcome of simulations throughout the fluctuation region. We also present all 
the data in Table I. The relative accuracy is high far from the fluctuation region (better then 1 %), and gets worse in 
the vicinity of the critical point where finite-size corrections are the largest. We show errorbars in all plots. 

First, we check for consistency between our approach and the previously obtained result for 9q, see Eq. I|2.17[) . 
Within the errorbars, the agreement is perfect, see Fig. 0] 

Knowing C and 8 is all we need to address the asymptotic MF behavior at large \X\, see Eqs. ()3.26|) . I|3.17jl . [|3.20JI . 
and 1)3. 23|) . The agreement with the MF theory based on the notion of quasicondensate is remarkable. Despite the 
fact that in MF we keep consistently only large terms cx X and oc \[X^ for numeric reasons the constant term also 
happens to be accurate. The self-consistent theory of Ref. |2J is claimed to go beyond conventional MF, but it is not 
known at present whether it reproduces correctly the vX terms. 

The agreement between the data and MF predictions makes it easy to estimate the size of the fluctuation region 
(in terms of X) to be about ~ 0.4 on the superfluid side of the transition, and roughly two times smaller, ~ 0.2, on 
the normal side. 

Other quantities of interest are the universal amplitudes for the superfluid and condensate densities in the critical 
regime, when f s = A S X", and /o = AoX" 1 - 1 ^. Here v = 0.6715 and r] — 0.038 are the correlation length and the 
correlation function critical exponents of the XY-universality class in d = 3, see, e.g., Ref. I2H By fitting the data for 
f s and /o at small X to the power laws with known exponents we obtain 

A s = 0.86(5) , (5.1) 



A = 0.89(5) . (5.2) 

^-effect. It appears in Fig. @ as if the total density has a cusp at the critical point. However, n is not singular at 
X = 0, and what looks as a cusp in Fig. (J2J) is in fact a relatively sharp crossover slightly shifted to the normal side 
of the transition point, see Fig. 

The crossover itself is predicted by MF equations since 0(X) changes its slope from roughly —X to X, see Eqs. Q3.2JI 
and (|3.26|) . We are not aware of any special reason why it has to be so sharp and so close to the transition point — the 
slope changes at X m —0.01 and the crossover region is only about AX <~ 0.01 in width. We call this surprising 
non-perturbative result the V-cffcct" as suggested by the shape of the ^-function. 



VI. DISCUSSION 



To discuss results in the canonical ensemble setup we need to change from the chemical potential as a control 
variable to reduced temperature. First, we note that our results immediately generalize to the case when the properly 
rescaled reduced density 

Y = 1 " c(r) (6 1) 

m 3 T 2 U n c (T) ' 1 ' ' 
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FIG. 4: The 8- and A-functions in the vicinity of the critical point, 
is used as a tuning variable. Indeed, 

Y = X(X) 
n c „d = m 3 T 2 Uf (X) 



(6.2) 



is nothing but the parametric dependence of the condensate density on density at a fixed temperature [and similarly 
for f s (Y) and 9(Y)}. 



- \ 



0.5 



1.5 



FIG. 5: The effective chemical potential, 9, dependence on reduced temperature variable Z. 
Next, we observe that Eq. (|2.9|l can be also considered as an equation for T c {n) 



3.313 



3/2 



Cm 3 T?U = n 



(6.3) 



Using it in the definition of the A-function we obtain 



mT c 
3.313 



3/2 



- Cm 3 T?U - 



mT 
3.313 



3/2 



Cm A T z U = m 6 T z U\{X) . 
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FIG. 6: The condensate and superfluid density dependencies on reduced temperature variable Z. 



Finally, keeping only the leading linear in U and t terms we arrive at 

Z = X(X) . (6.4) 

where Z is the rescaled reduced temperature parameter 

„ 3n 



2m 3 T 2 U 



t . (6.5) 



This establishes the functional relation between t and X, and the parametric dependence of other properties on t. 

In Figs. ISJand we plot 9, fo and f s dependencies on the reduced temperature variable Z. From Fig. [5]we deduce 
the size of the fluctuation region consistent with the previous estimate in terms of X and relation Z = X(X), i.e. it 
is about ~ 0.5 on the superfluid side and ~ 0.1 on the normal side. Amazingly, in temperature plots for n cn d(r) and 
n s (T) one may not even distinguish between the MF theory and numerical data on large scale; only in the inset which 
covers just the fluctuation region we can see that the MF curve goes wrong very close to T c . 

We also verify the prediction of the self-consistent theory 2 ^ for the critical amplitude of the condensate density. 
As mentioned previously, the self-consistent theory makes predictions for the critical region as well; the relation 
between the condensate density and reduced temperature is derived as [see Eq. (18) in Ref. 0] fo ~ 0.251 Z v 
with v 1 = 1/(1 + yJh/Qii) « 0.66. The condensate density exponent is accurate, though the derivation was done 
assuming that the self-energy behaves as S(fc — > 0) ~ fc 3 / 2 instead of the true critical behavior — > 0) ~ A; 2-17 with 
r) = 0.0380(4) (see Ref. l25l and the discussion below). Since d\/dX\x=o — \' ~ 1-8(1) the amplitude in the critical 
law f = A T) is directly related to the value of A a in Eq. (JOJ 

A { P = Ao ■ « 0.62(4) , (6.6) 

in disagreement with the theory by a factor of two. 

To characterize the fluctuation region in terms of the gas parameter n 1 / 3 a we identically rewrite 

Z = 91.95nV3 a (6J) 

The coefficient in the denominator is too large to be ignored even in qualitative estimates. It was mentioned in 
Ref. that Bose gases should demonstrate universal properties only at na 3 <C 10~ 5 . This statement is unambiguously 
confirmed by the present study, since for larger values of na 3 the fluctuation region in temperature is already of order 
T c itself. The other identical way of writing Z is 

0.249 k T , . 

Z = — -t . (6.8) 
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Obviously, the idea of universality is meaningful only if it is possible to separate ideal-gas short-wavelength physics from 
strongly coupled long- wavelength fluctuations, i.e. when k c -C fcy, and thus, according to Eq. I|6.8[) . the fluctuation 
region in reduced temperature is small. 

In a typical experiment with ultracold gases the parameter na 3 is as small as ~ 10~ 6 , and the system may be 
considered as weakly interacting. The interaction induced critical temperature shift in the homogeneous system is 
very small (the corresponding critical density shift is given by Eq. I|2.9|) and is only about 1% for fc c /fcr ~ 1) and is 
difficult to study experimentally. This does not imply, however, that all non-perturbative effects in the vicinity of the 
critical point are of academic interest as well. 

The quantity directly relevant to the experimental setup is X(X) since it describes, according to Eq. i|2.5|l . the 
density profile in the trapping potential if it is smooth enough to guarantee the hydrostatic regime. In this regime 
the density variation over the mode-coupling radius r c ~ l/fc c reduces to n = n{T 1 /i(R)), where n(T, /i) is the 
homogeneous equation of state, /i(R) = const — V^ x t(R), and V ex t is the trapping potential. It follows from Figs. |3 
and 0] that the "cusp" on the density profile does not coincide with the point where the condensate first appears, but 
is slightly shifted to the perimeter of the trap. Moreover, by fitting experimental data away from the "cusp" using 
MF equations, one may directly measure the interaction induced universal chemical potential shift [through 9q + 2C\. 

The variation of the A(X)-function across the fluctuation region is about ~ 0.3, which transforms into density 
variation of order 

An/n c ~ 40n 1/3 a , (6.9) 

or a 40 % strong effect for na 3 = 10~ 6 ! In Fig. [7] we show the density profile in a smooth parabolic trap 

V ext = mtu 2 R 2 /2 = T(R/R T ) 2 , (6.10) 

when the condensate density in the middle of the trap is comparable to the critical density. Amazingly, in this 
situation the fluctuation region extends all the way from the critical point to the trap center, and the MF theory 
completely fails to describe the superfluid side. Thus, the non-perturbative physics of the fluctuation region and the 
prediction that it is completely characterized by the classical field theory can be studied even by experimenting with 
very dilute gases, na 3 < 10 -6 . As Fig. [7| clearly demonstrates, all effects in the middle of the trap are strong. 

It is also worth mentioning that in the fluctuation region the density profile can not be decomposed into the sum 
of the smoothly varying, monotonic non-condensate density background and the condensate density bump. The 
normalized condensate density in Fig. [7| increases faster then n/n c , and the naive "decomposition" technique would 
underestimate the actual condensate density by almost a factor of two. 

For the quasi-homogenenous description to work, it is necessary to keep the external potential gradients small. If 
the chemical potential in the middle of the trap corresponds to (fx — /i c ) = m 3 T 2 U 2 X(0), then the critical point is 
located at a distance R c = (mTU/uS) \j2Xifij. The hydrostatic approach can be used when k c R c ^ 1, or 

~ < m 3 TU 2 ^2X{0) sa 740(na 3 ) 2/ V^(0) • (6.11) 

For the parameters in Fig. 0we need then u <C 0.04T — a condition easily satisfied experimentally 26 . The crucial 
point is then in achieving the necessary spatial resolution in shallow traps. 

Our final result concerns the universal part of the self-energy at the critical pointjiiS* 2 *. The most recent calculation 
based on renormalization group equations for the vertex functions^ predicted {x = k/k c ): 

a{x) = ^ [E(fc c x) - E(0)] ^ Sx 2 ^ , (6.12) 

with S = 1.54, n w 0.104, and an extended crossover between the free particle and critical regimes (we preserve the 
present paper definition of k c — m TU). 

This quantity is directly related to the occupation numbers in the k — ■> limit by (k 2 /2mT) = 1/[1 + a(x)/x 2 ]. 
In Worm algorithms, the one-particle density matrix is calculated automatically as the central part of the numerical 
schemsii, and thus — ^ r ,o(r)e lkr is readily available. In Fig. [S] we plot the universal part of the occupation 
number distribution as ln[2mT/(fc 2 nk)} versus \n(x) (in the inset we plot ln[T/(ek nk)] versus ln(x), i.e. subtracting 
the leading bare spectrum dependence) for the system size 128 3 and see the crossover between the free particle, 
1/rik (x x 2 , and the U(l) critical, cx a; 2-17 , behavior at x ~ 1. The smallest values of x are effected by finite size 
effects and axe not shown in the inset (in finite systems ti^—q — Hoid rs " > 

1/L 1+ '' is finite at the critical point). 
The correlation function exponent is known very accurately 2 ^ r\ « 0.0380(4), and we consider it as known (our 
data are consistent with this value). Clearly, very little changes occur in the structure of the k 2 /2m + S(fc) — £(Q) 
expression across the fluctuation region. It seems the best way to characterize the crossover is to write the whole 
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FIG. 7: The normalized density profile of the weakly-interacting atomic gas with na 3 — 1CP 6 in a smooth harmonic trap 
(circles). The chemical potential in the middle of the trap corresponds to X(R = 0) = 0.3 (see text). The dashed line is the 
mean-field theory prediction. We also show the normalized condensate density (squares) for comparison. 




FIG. 8: The universal part of the occupation number distribution ln[2mT / 'fe^rik] at the critical point. Scattering of points is 
due to statistical errorbars. The dashed line in the main plot is the bare dispersion law contribution, ln[e k /(fc;?/2m)], and the 
solid line is the critical law (2 — r/)x + \nS. In the inset we plot ln[T/ekrok] to see the crossover more clearly. Note the change 
in the vertical scale. 

expression as x 2 + a(x) = x 2 ^^ x ^ with f{x) interpolating between and rj = 0.038. Also, the crossover is localized 
in the vicinity of x ~ 1. The numerical value of S in the asymptotic limit is found to be very close to unity 

S = 1.038(6) (6.13) 
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TABLE I: Final results for the universal functions 0{X), X(X), fo(X), and f s (X). 


X 


6{X) 


\{X) 


MX) 


fs(X) 


-3.738 


3.17 


-0.406(1) 






-3.338 


2.82 


-0.381 






-2.938 


2.47 


-0.355 






-2.538 


2.13 


-0.327 






-2.138 


1.79 


-0.297(1) 






-1.738 


1.45 


-0.263 






-1.498 


1.26 


-0.242 






-1.338 


1.13 


-0.226 






-1.178 


1.00 


-0.210(1) 






-1.018 


0.874 


-0.193 






-0.8584 


0.751 


-0.174 






-0.5384 


0.517 


-0.130(1) 






-0.3064 


0.361 


-0.0875(3) 






-0.2424 


0.323 


-0.0804(2) 






-0.2104 


0.304 


-0.0732(4) 






-0.1544 


0.275 


-0.0593(6) 






-0.1384 


0.267 


-0.0546(7) 






-0.1304 


0.263 


-0.0527(2) 






-0.1064 


0.252 


-0.0462(3) 






-0.09039 


0.246 


-0.0409(3) 






-0.06639 


0.237 


-0.0324(5) 






-0.05839 


0.234 


-0.0299(6) 






-0.04479 


0.229 


-0.0251(7) 






-0.03199 


0.227 


-0.0213(12) 






-0.01839 


0.225 


-0.0173(11) 






-0.01039 


0.226 


-0.0134(7) 






-0.002393 


0.235(1) 


-0.0046(18) 






0.001607 


0.244(2) 


0.0021(15) 


0.0116(30) 


0.0130(30) 


0.005607 


0.255(2) 


0.0094(18) 


0.0240(15) 


0.0259(16) 


0.009607 


0.266(2) 


0.0171(11) 


0.0356(10) 


0.0381(11) 


0.01361 


0.277(1) 


0.0243(8) 


0.0446(7) 


0.0482(8) 


0.02161 


0.294(3) 


0.0367(23) 


0.0630(6) 


0.0674(7) 


0.02850 


0.307(2) 


0.0456(8) 


0.0772(8) 


0.0826(8) 


0.04350 


0.332(2) 


0.0663(8) 


0.105(1) 


0.111(1) 


0.05361 


0.349(1) 


0.0802(21) 


0.123 


0.132 


0.07761 


0.394(2) 


0.117(3) 


0.160(2) 


0.173(2) 


0.08561 


0.406(1) 


0.127(2) 


0.173(1) 


0.187(2) 


0.1016 


0.431(1) 


0.147(4) 


0.197(2) 


0.215(1) 


0.1256 


0.469(1) 


0.178(3) 


0.237(4) 


0.255(3) 


0.1416 


0.492(2) 


0.197(4) 


0.260(2) 


0.279(3) 


0.1576 


0.518(2) 


0.219(3) 


0.287(4) 


0.310(6) 


0.1816 


0.555(1) 


0.249(2) 


0.320(3) 


0.342(4) 


0.2089 


0.587(8) 


0.276(4) 


0.360(8) 


0.386(5) 


0.2376 


0.637(3) 


0.317(4) 


0.401(4) 


0.426(9) 


0.2616 


0.672(3) 


0.348(5) 


0.434(5) 


0.462(4) 


0.2936 


0.716(2) 


0.385(3) 


0.479(2) 


0.511(6) 


0.4216 


0.895(2) 


0.539(3) 


0.652(3) 


0.683(7) 


U.DDlD 


1.21(1) 


u.oio^y ) 


u.yoD^ / j 


U.yyo^lU ) 


1.062 


1.72 


1.27(1) 


1.44(1) 


1.51(1) 


1.462 


2.22(1) 


1.72(1) 


1.93(1) 


2.00(2) 


1.862 


2.69(1) 


2.15(2) 


2.40(1) 


2.47(3) 


2.662 


3.63(2) 


3.01(3) 


3.28(4) 


3.39(7) 


3.462 


4.56(3) 


3.89(4) 


4.20(3) 


4.33(5) 
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